clear all; close all; clc;

% Внутренние параметры
R = 120; % Ом
C = 33e-12; % Ф
L = 100e-6; % Гн
C1 = C/2;
R1 = 2*R;
Rd = 0; % Параметры диода
rd = 0;
Uo = 0;

omega0 = 1/sqrt(L*C1); % Резонансная частота
f0 = omega0 / 2 / pi;
fprintf('f0 = %f MHz\n', f0/1e6);

% Ось времени
T = 1 / (10000*f0);
t = 0:T:(15 * 1/f0);
lt = length(t);

% Выделение памяти
Uc = nan(1, lt);
Ul = nan(1, lt);
dUc = nan(1, lt);
Ud = nan(1, lt);
i = nan(1, lt);
di = nan(1, lt);

% Внешние параметры, внешние фазовые переменные
E = 10*ones(1, lt);

% Начальное состояние
Uc(1) = 0; % Начальное напряжение конденсатора
i(1) = 0; % В начальный момент ток в цепи отсутствует
di(1) = 0;
Ul(1) = di(1) * L; % Напряжение на катушке

for k = 2:lt
    % Напряжение на диоде
    if (E(k-1) - Ul(k-1) - Uc(k-1)) < Uo
        Ud(k-1) = Rd*i(k-1);
    else
        Ud(k-1) = rd*i(k-1) + Uo;
    end
    
    Ul(k) = E(k) - Uc(k-1) - i(k-1)*R1 - Ud(k-1); % Напряжение на катушке
    di(k) = Ul(k) / L; % вызывает приращение тока,
    i(k) = i(k-1) + di(k-1)*T;
    dUc(k) = i(k) / C1;% который определяет заряд конденсатора
    Uc(k) = Uc(k-1) + dUc(k)*T;
end

U = Uc / 2;

figure(1);
plot(t*1e9, [Uc; E; U]);
xlabel('t, ns')
ylabel('U_c, E, U, Volt')
legend('U_c(t)', 'E(t)', 'U(t)')
grid on



